# comment to keep archive current

from pycuda.compiler import SourceModule


MatchTextureModule = SourceModule("""

#include <float.h>

texture<float, 2, cudaReadModeElementType> ref_texture;
texture<float, 2, cudaReadModeElementType> img_texture;
texture<float, 2, cudaReadModeElementType> wt_texture;


__global__ void CopyTexture(float * image, int N)
	{
	int rshift = N*N*blockIdx.x*gridDim.y;

	for (int i=0; i<N; i++) {
		int rshift1 = i*N*gridDim.y + N*blockIdx.y;
		int row = N*blockIdx.x + i;
		for (int j=0; j<N; j++) {
			int col = N*blockIdx.y + j;
			image[rshift + rshift1 + j] = tex2D(ref_texture,col,row);
			}
		}
    }


__global__ void CopyTexture1(float * image, int N, int P)
	{
	for (int i=0; i<N; i++) {
		for (int j=0; j<P; j++) {
			image[i*P + j] = tex2D(ref_texture,j + 0.5f,i + 0.5f);
			}
		}
    }


__global__ void ShiftTexture(float * image, int N, float sx, float sy)
	{
	int startx = max(0,int(floor(sx)));
	int stopx = min(N,startx+N);

	int starty = max(0,int(floor(sy)));
	int stopy = min(N,starty+N);

	for (int i=startx; i<stopx; i++) {
		for (int j=starty; j<stopy; j++) {
			image[N*i + j] = tex2D(img_texture,i-startx,j-starty);
			}
		}
	}

__global__ void ProjectionMatch(int N, int K, int * best_idx, float * resid)
	{
	float best_resid=FLT_MAX;
	int best_index=0;

	for (int k=0; k<K; k++) {
		int2 begin, end;
		float the_sum=0.0f;
		float3 wtran;
		wtran.x = tex2D(wt_texture,0,k);
		wtran.y = tex2D(wt_texture,1,k);
		wtran.z = tex2D(wt_texture,2,k);

		begin.x = -int(floorf(wtran.x));
		end.x = min(begin.x + N,N);
		end.x -= int(ceilf(wtran.x)) - int(wtran.x);
		begin.x = max(0,begin.x);

		begin.y = -int(floorf(wtran.y));
		end.y = min(begin.y + N,N);
		end.y -= int(ceilf(wtran.y)) - int(wtran.y);
		begin.y = max(0,begin.y);

		wtran.x += 0.5;
		wtran.y += 0.5;

		for (int i=begin.x; i<end.x; i++) {
			float xpos = blockIdx.x*N + i + wtran.x;

			for (int j=begin.y; j<end.y; j++) {
				float ypos = blockIdx.y*N + j + wtran.y;
				float diff = tex2D(ref_texture,ypos,xpos)
								- tex2D(img_texture,N*threadIdx.x + j,i);
				the_sum += diff*diff;
				}
			}


//		the_sum /= wtran.z;


		float swapf = best_resid - the_sum;
		swapf = fminf(1.0f,swapf);
		swapf = fmaxf(0.0f,swapf);
		swapf = ceilf(swapf);
		int swapi = int(swapf);

		best_resid = (1.0f - swapf)*best_resid + swapf*the_sum;
		best_index = (1 - swapi)*best_index + swapi*k;
		}

	int idx = threadIdx.x*gridDim.x*gridDim.y;
	idx += blockIdx.x*gridDim.y + blockIdx.y;
	resid[idx] = best_resid;
	best_idx[idx] = best_index;
	}



__global__ void ProjectionMatch1(int N, float * resid)
	{
	float the_sum = 0.0f;
	for (int i=0; i<N; i++) {
		float row = N*blockIdx.x + float(i);
		for (int j=0; j<N; j++) {
			float col = N*blockIdx.y + float(j);
			float diff = tex2D(ref_texture,col,row)
							- tex2D(img_texture,N*threadIdx.x + j,i);
			the_sum += diff*diff;
			}
		}

	int idx = threadIdx.x*gridDim.x*gridDim.y;
	idx += blockIdx.x*gridDim.y;
	
	resid[idx + blockIdx.y] = the_sum;
	}
	""")

